#!/usr/bin/octave -qf
sizes = [35, 50, 75, 100, 139, 150, 160, 200, 250, 300, 340, 350, 400, 450, 490, 500];

arg = argv;
sample = load("-ascii", arg{1});

nbins = 200;
bins = min(sample)-1:floor((max(sample)-min(sample))/nbins):max(sample);
hist = histc(sample, bins);
ix = find(hist < 50 & hist > 2); # the numbers depend or nbins
snoise = bins(ix(1)); # everything that occurs more frequently must be noise
sample(find(sample < snoise)) = 0;

order = 2;
signal = dwfilt(sample, order); # Savitsky-Golay with Durbin-Watson criterion, order = 2
d1 = dwfilt(diff(signal), 1);
d2 = dwfilt(diff(d1), order);
d3 = dwfilt(diff(d2), order);

signal = signal(3:length(d3)+2);
d1 = d1(2:length(d3)+1);
d2 = d2(2:length(d3)+1);

thr00 = snoise;
thr0 = snoise;
thr1 = 3*noise(d1);
thr2 = 3*noise(d2);

region = find(d2 < -thr2 & signal > thr0 & (d1 > thr1 | d1 < -thr1) );
mark = zeros(1, length(signal));
mark(region) = ones(1, length(region));

% bridge the gaps (if two marked regions are one point apart
% it's likely because the d1 constraint above knocked out a point near zero)
for i = 2:length(mark)-1
  if (mark(i) == 0 & mark(i-1) == 1 & mark(i+1) == 1)
    mark(i) = 1;
  endif
endfor

starts = find(diff(mark) > 0);
ends = find(diff(mark) < 0);
common_size = min(length(starts), length(ends));
starts = starts(1:common_size);
ends = ends(1:common_size);

% extend the peak regions by annexing the areas of substantial 1st derivative
region = {}; # reuse the region array for a different purpose
for i = 1:length(starts)
  from = starts(i);
  to = ends(i);
  if to - from <= 3 # discard short regions as noise
    continue;
  endif

  fusedbefore = fusedafter = 0;

  j = from;
  while ( d1(j) > thr1 & signal(j) > thr0 )
    if ( !fusedbefore & sign(d3(j)) != sign(d3(j-1)) )
      fusedbefore = j;
    endif
    j--;
  endwhile
  from = j;
  % don't bother with small fused areas
  if fusedbefore - from <= 7
    fusedbefore = 0;
  endif

  j = to;
  while ( d1(j) < -thr1 & signal(j) > thr0 )
    if ( !fusedafter & sign(d3(j)) != sign(d3(j-1)) )
      fusedafter = j;
    endif
    j++;
  endwhile
  to = j-1;
  if to - fusedafter <= 7
    fusedafter = 0;
  endif

  if ( fusedbefore )
    region{end+1} = [from, fusedbefore-1];
    from = fusedbefore;
  endif
  if ( fusedafter )
    region{end+1} = [from, fusedafter];
    region{end+1} = [fusedafter+1, to];
  else
    if to - from > 5
      region{end+1} = [from, to];
   endif
 endif
endfor

% this detection routine somehow resuls in overlapping
% regions. Join those that overlap by more than one point
region1 = region;
region = {};
i = 1;
while ( i < length(region1) )
  reg1 = region1{i}(1):region1{i}(2);
  reg2 = region1{i+1}(1):region1{i+1}(2);
  newreg = intersect(reg1, reg2);
  if length(newreg) >= 2
    newreg = union(reg1, reg2);
    region{end+1} = [newreg(1), newreg(end)];
    i = i + 2;
  else
    region{end+1} = [reg1(1), reg1(end)];
    i++;
  endif
endwhile

if length(intersect(region{end}(1):region{end}(2), region1{end}(1):region1{end}(2))) <= 1
  region{end+1} = region1{end};
endif

for i = 1:length(region)
  range = region{i};
  from = range(1);
  to = range(2);
  [mind2, mind2x] = min(d2(from:to));
  sizex = from + mind2x - 1;
  if mind2 > -thr2 # this is not a real peak
    sizex = floor((from + to)/2);
  else
    % stop toward real maximum if there is an offest in the 2nd derivative
    if sizex != from & sizex != to
      if signal(sizex - 1) > signal(sizex)
	sizex = sizex - 1;
      endif
      if signal(sizex + 1) > signal(sizex)
	sizex = sizex + 1;
      endif
    endif
  endif

  edge1 = signal(from)/2;
  edge2 = signal(to)/2;
  base = edge1 + (sizex - from)*(edge2 - edge1)/(to - from);
  height = signal(sizex) - base;
  if height < thr00 # height small or negative? go away.
    continue;
  endif
  basearea = mean([edge1, edge2])*(to - from);
  area(i) = sum(signal(from:to)) - basearea;
endfor

height = zeros(1,length(starts));
pos = zeros(1,length(starts));
for i = 1:length(starts)
  from = starts(i);
  to = ends(i);
  [height(i), p] = max(signal(from:to));
  pos(i) = from + p;
endfor

n = length(sizes);
% take the last n peaks and check that they their height is similar
pos = pos(end-n-1:end);
height = height(end-n-1:end);
m = mean(height)
sd = std(height)

for i = 1:length(pos)
  printf("%g\t%g\t%g\n", pos(i), height(i), area(i));
endfor
